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ABSTRACT 

The extreme eigenvalues of the Chebyshev pseudospectral differentiation operator are 
0 (A 2 ) where A is the number of grid points [4]. As a result of this, the allowable time 
step in an explicit time marching algorithm is 0 (A -2 ) which, in many cases, is much below 
the time step dictated by the physics of the P.D.E. In this paper we introduce a new set of 
interpolating points such that the eigenvalues of the differentiation operator are 0 (A) and 
the allowable time step is 0 (A -1 ). The properties of the new algorithm are similar to those 

of the Fourier method but in addition it provides highly accurate solution for nonperiodic 

! 

3 

boundary value problems. 
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1. Introduction 


Consider the first order hyperbolic initial boundary value problem 

u t — u x = 0 — 1<x<1,£>0 (1.1) 

u(x,0) = u°(x) — 1 < x < 1 (1.2) 

u(l,t) = s(t ) t > 0. (1.3) 

A standard pseudospectral method [5] for solving (1.1)-(1.3) is based on interpolating the 
solution at the extremal points 

*< = cos (£) * = 0, • • • , AT (1.4) 

of the Nth-order Chebyshev polynomial 

Tjv(x) — cos(N arccos(x)). (1.5) 

Using this method for space discretization and a standard explicit scheme (e.g., Runge- 
Kutta) for time discretization, one encounters a stability condition which has to satisfy [5] 

At = 0(N- 2 ). (1.6) 


This restriction is very stringent and forces the user to march in time steps which, in many 
cases, are much bellow the time step dictated by the physics of the problem. A way of 
overcoming this annoying phenomenon is to use implicit or semi-implicit time marching 
techniques. Since the pseudospectral differentiation matrix is dense, the resulting algorithm 
is highly time consuming. Therefore, we would like to find a way to advance explicitly in 
time with a less restrictive stability condition. Our research is aimed at this target. 

Chebyshev points (1.4) are bunched near the boundaries with minimal spacing of 0(N~ 2 ). 
Since the pseudospectral method is global, there is no direct relation between the minimal 
spacing and the stability condition as in the finite-difference method [9]. Nevertheless, nu- 
merical experience and heuristic reasoning led us to ‘blame’ the super fine grid near the 
boundaries for the severe stability condition (1.6). When there are sharp gradients near the 
boundaries, the clustering of points is needed for resolution and the small time step can 
also be anticipated by physical reasonings. But when the high gradients are elsewhere, or 
if the solution is evenly smooth, there seem to be no justification in putting more points 
near the boundary. Thus, we are led to the conclusion that the numerical tool we are using 
(polynomial interpolation) is not appropriate in these cases. 

In [2], Bayliss et al. describe a physical model with very sharp gradients. In order to 
overcome the numerical difficulties, they have designed an algorithm where the problem is 


1 



transformed so as to minimize some functional. In our research, we use a similar transfor- 
mation approach. The collocation points are chosen as 

«) - 1 < < 1 » = o < a < l, (i. 7 ) 

where 

= cos^j i = (1.8) 

a is a parameter and g(y, a) is a ‘stretching function’. By a proper choice of the parameter 
q we increase the minimal spacing near the boundaries such that 

Aar min = 0 (N- 1 ) . (1.9) 

Consequently, we are able to advance in time with the favorable stability condition 

At = 0(N~ 1 ). (1.10) 

Moreover, it is shown that, as N — + oo , one needs only two points per wavelength for 
resolution (as in the Fourier method ) and not 7r points as in the Chebyshev case [5]. Thus, 
fewer points are needed to model the P.D.E. (a saving of almost 40% ). The transformation 
function is described in Section 2. In Section 3 we present the resolution analysis which also 
reveals the approximation subspace on which the solution is projected. In order to efficiently 
implement the algorithm, it is important to choose the appropriate parameter a and this 
subject is discussed in Section 4. In Section 5 we present a more general transformation 
which gives additional flexibility to the new interpolation method. The paper is concluded 
in Section 6 in which we present numerical results. 

2. New Interpolation Points 

Chebyshev pseudospectral solution of (1.1)-(1.3) is based on approximating the spatial 
derivative by differentiating analytically the interpolating polynomial. If v is an N 
dimensional vector which approximates u(x) at the interpolation points (1.4) then the 
vector 

v = Dv (2.1) 

approximates u'(x) at (1.4). D is the spatial differentiation matrix which incorporates the 
boundary condition (1.3). The entries of D are given in [6] ((2.1) can be accomplished by 
using FFTs requiring only O(NlogN) operations [6]). The matrix D is very ill-conditioned, 
with eigenvalues scattered in the left side of the complex plane [4]. While most of the 
eigenvalues grow like O(N), a few of them are 0(N 2 ) [4]. These extreme eigenvalues 
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can be considered as the reason for the severe stability condition (1.6). We have to choose 
At — 0(N~ 2 ) so that all the eigenvalues of A tD will be included in the domain of stability 
of the time marching scheme. 

Furthermore, 


Ax min = min |x,+l - x,| = 1 - cos(7 t/N) - 0{N 2 ). 

t 


( 2 . 2 ) 


This phenomenon of having domain of eigenvalues whose size is proportional to the reciprocal 
of the minimal spacing is typical to many differentiation matrices. Even in cases where this 
correspondence does not hold, as in the Legendre pseudospectral method [10] , we still have 
to choose time step dictated by the minimal spacing due to numerical instability whose origin 
is the ill-conditioning of the matrix which diagonalizes the differentiation matrix [12]. Thus, 
we would like to have a set of interpolating points with larger minimal spacing. We are going 
to attain this goal by mapping the Chebyshev points (1.8) to another set of points in [—1, 1] 
such that the minimal spacing near the boundary is ‘stretched’. 

Let us consider the transformation 


= g(y,a) = 


arcsin(ay) 


arcsin(a) 

Computing the derivative at the grid points X{ 


x,V € [-!,!]■ 


(2.3) 


Xi = g(yi\a) 0 <i<N, 


(2.4) 


where 


y,- = cos ^ 0 < i < N (2-5) 

is accomplished by making use of the chain rule. For given / € C l [— 1, 1] , we have 


Hence, we modify (2.1) to read 
where A is a diagonal matrix 


d£ _ 1 d/ 

dx g'(y;a)dy' 

v — ADv 
1 


An — , 


g'(yi\ a ) 


and 


g'(y\ a ) = 


a 


( 2 . 6 ) 

(2.7) 

( 2 . 8 ) 

(2.9) 


arcsin(a) ^1 _ (ay) 2 ’ 

v and v' contains the approximated values of u(x) and u (x) respectively at x; = g{yi\a). 
We have 
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Lemma 2.1. If x,- , 0 < i < N, satisfy (2-4) then the minimal spacing between the points 
is attained near the boundaries. 


Proof: Define 


Using the mean value theorem 


0 = arccos(y). 


( 2 . 10 ) 


Ax* — i — 

= ^(6) M, 0i<Zi<0i+ 1, 0 < i < TV — 1 


while 


By (2.3) and (2.10) we have 


where 


Since 


dg 


Ad = — . 

N 


a 


h{6) = 


dd arcsin(a) 
sin($) 


m 


\J\ — (acos(0)) 3 


( 2 . 11 ) 

( 2 . 12 ) 

(2.13) 

(2.14) 


h'(0) = cos ffl( a2 ~ 1) 


W = 7 — S’ (2.15) 

(l - (acos(0)) 3 )’ 

is negative in 0 < d < | and positive in | < d < tt , — h(d ) attains its minima at d = 0,7r 
and the result follows. 

Furthermore, 


Ax r 


= l-x l arcsin (<* <*«(£)) 


• , , ■ (2.16) 

arcsin(o:) v ; 

It is easily verified that the R.H.S. of (2.16) is monotonically increasing with a , a € (0, 1) 
and 

lim Ax m j n = — (as in Fourier case) (2-17) 

„iJo ^ Xmin = ^ — cos ( as i n Chebyshev case). (2.18) 


Lemma 2.2. If 


then 


a = 1 -Jn +0 ( N ~ 3 ) 


a 2tt 

— /— ~ 
\Ar 2 + 2c + 


wM +0 ^ ■ 


(2.19) 

( 2 . 20 ) 
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Proof: Define z = jj then 


where 


AXmin = 1 ~ h(z ) 

_ arcsinfofc)) 

• ( ( \\ ? 

arcsin($r 2 (z)) 

flTi(ar) = (1 - cz 2 ) cos(7rz) + 0 (z 3 ) , 
52(2) = (1 - cz 2 ) + O (z 3 ) . 


We have h( 0) = 1 and 


k'(0) = - lim 


ffi 


9 2 


9\ t 


Using L’Hospital’s rule it easily verified that 

h\ 0) = -;(Vft(0)-VS(0)) 

= ^ (Vir 2 + 2c - >/&) . 


( 2 . 21 ) 

( 2 . 22 ) 

(2.23) 

(2.24) 

(2.25) 


(2.26) 


and the result follows. Based on Lemma 2.2 we conjecture that the time restriction of the 
new interpolation method satisfies (1.10). Numerical results reported in section 6 assist this 
conjecture. 

3. Resolution Analysis 

Let 

fi(x) = cos(r7rx) — 1 < x < 1, (3-1) 

/ 2 (x) = sin(r7rx) — 1 ^ 2: < 1 (3-2) 

be functions whose derivative we want to approximate ( r is a real number which indicates 

the wave number). Substituting (2.3) in (3.1) and (3.2) we have 

fi(x) = h(y) = cos[ro arcsin(ay)] (3.3) 

h(x) = f 2 (y) = sin[marcsin(c«/)] (3.4) 

where 


V7T 


m = 


arcsin 


in(a) 


(3.5) 


In the appendix we show that 
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for m even 

cos(m^) = (-l)“7 1 m (sin(V»)), 

(3.6) 

for m odd 

sin(m^) = (-l)VT m (sin(0)). 

(3.7) 

Hence, using (3.3)-(3.5) 

fi(y) = (-l)“T m (aj/) ( m even ) 

(3.8) 

and 

My) = (-1) 125-1 T ro (ay) ( m odd). 

(3.9) 

T m (ay) is a polynomial 

in y, therefore, interpolating at N + 1 

points, m < N, will 


result in the function itself. Hence, the new algorithm is exact for the following N + 1 
functions (N even) 

1, cos(2px), cos(4px), ■ • • , cos(Npx), sin(px), sin(3px), • • ■ , sin ((N — l)px) (3.10) 

where 

p = arcsin(a). (3.11) 

The set of functions (3.10) span the approximation subspace. An elaborate discussion of this 
subspace will be given in future paper. 

Let us clarify now what we mean by resolution. If the Chebyshev expansion of a function 
h{y) is 

0° J 

M y) = X] ~ a kTk(y) Co = 2, c k = 1 for k ^ 0, (3.12) 

k=0 Ck 

and there is K such that a k decreases rapidly when k increases beyond K then we say 
that h(y) is resolved by K terms. Since Chebyshev expansion is qualitatively similar to 
interpolation in Chebyshev points (1.8), it is equivalent to saying that we need K points 
in order to resolve h(y) by interpolation. 

We speculate that 

Conjecture: The function T m (ay ) is resolved by M + 1 terms where 

M = [am] . (3.13) 

The reasoning for this conjecture goes as follows: let 

m i 

T m (ay) = £ -a?T k (y) cq = 2, c k = 1 for k ± 0, (3.14) 

k = o Cfc 
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then 


(3.15) 


m _ 2 r T m ( a y)T k (y) 

H ~ * l v'T^F 

Chebyshev polynomials satisfy the recurrence relation 

W*) = 2 xTn(x) -r n _!(x). (3.16) 

Hence 

a m = l f { 2 ayr ro -i(ay) - T m - 2 (ay)} r t (y) ^ m > 2 ( 3 . 17 ) 

k wJ yfi - y 2 

a° = 0 , aj = 0 , a} = a. (3.18) 

By (3.16), 2yT k {y) = T fe+ 1 (y) + T k -i(y) which imply the relation 

a? = c k a (a^Ji 1 + a^ 1 ) - <*' 2 k > 0 , m > 2, (a ™ -1 = 0). (3.19) 

We have programmed (3.19) with initial values (3.18) and ran it for many values of m and 
a and have always observed that when k< M , af are nondecreasing. Once k is greater 
then M , a™ decreases very rapidly. 

Therefore, using (3.5), the maximal wave number which can be resolved by the new 


method is 

N arcsin(a) 

^ max — 

7ra 

(3.20) 

while, since (3.10) 

N arcsin(a) 

r N - 

7T 

( 3 . 21 ) 

is the largest wave 

number of mode resolved exactly by the new method. 


Observe that 

lim r max = N/2 (as in Fourier case), 

(3.22) 


lim r max = N/n (as in Chebyshev case). 

cr — 0 

(3.23) 

Thus, by (3.22), asymptotically, two points per wavelength are needed for resolution. 

We have shown that for r satisfying (3.5) with m even(odd), fi(y) (/ 2 (y)) is a 
polynomial in y. Let us discuss now the resolution of general trigonometric function 


f(x) = exp(tnrx). 

(3.24) 

We have 

f*7T 

f(v) = f(9(y ;«)) - ex Pt l arcsin ( a ) arcsin («y)l- 

(3.25) 
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If 


^ ( - a ) = k + P . 0</?<l ,* integer (3.26) 

then 

/(y) = exp[ifcarcsin(cn/)] exp[j/?arcsin(ay)]. (3.27) 

Let us assume, without loss of generality, that k is even then by Lemmas A.l, A. 2 

f(y) = (~1) ? (Tkiay ) - i\Jl - {ay) 2 Pk-i(cty)^ (yjl - (ay) 2 + iayj . (3.28) 

Resolution of /(y) by interpolation is influenced by the degree k of the polynomials 
involved and by the singularities at ±£. The degree of the interpolating polynomial has to 
be at least ak but the asymptotic rate of convergence depends only on the singularities. 
The relevant theory is presented below in greater generality. 

Let f(x) be an analytic function in £7 Z> [—1,1]. Since y'(y; a) (2.9) has singularities 
at y = ±1 /a , so does /(y) = /(y(y;o)). Define 

B - H ■ i) ^ 

and assume that a is close enough to 1 such that 


Bcg~ l (E). (3.30) 

Thus, /(y) is analytic in B. The rate of convergence of polynomial interpolation at 
Chebyshev points is based on the following theory from [13]: let K be a bounded continuum 
in C such that K c = the complement K is simply connected in the extended plane and 
contains the point at infinity. For such K there exist a conformal mapping ^(u;) which 
maps the complement of the unit disc onto K c [13]. Let $(.z) be the inverse of ^(in) and 


B t = {z: |$(*)| = f} ( t > 1) 


(3.31) 


denote the level curves in K c . 


Theorem 2.1: Suppose t > 1 is the largest number such that f(z) is analytic inside B t . The 
interpolating polynomials P n {z) with interpolating points z" that are uniformly distributed 
on K then satisfy 


lim max \ f(z 

n-*oo z£K 1 v 



(3.32) 
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Since Chebyshev points satisfy the definition of uniformly distributed points on [-1,1] 
[13], the asymptotic rate of convergence can be computed by making use of this theorem. 
We choose K = [-1,1], and the relevant conformal mapping is given by [8] 

3>(t /) = y ± \fy 2 — 1 ■ (3.33) 


Assume now that if f(y) has additional singularities in the complex plane then a is close 
enough to 1 so that the largest t corresponds to the singular points ±^. Thus 

t _ 1 + aA ~ <* 2 (3.34) 

a 


and the asymptotic rate of convergence is 

1 _ 1 — y/l -a 2 
t a 

Hence, by interpolating at 1V + 1 points, the asymptotic accuracy is ce where 


(3.35) 


(3.36) 


and c is constant which depends on / but does not depend on N or y. 


4. On the Choice of the Parameter a 

For a predetermined degree N , we would like to choose the appropriate parameter 
a. We give below three constructive ways for choosing a, based on different considerations. 

Resolution considerations - Sometimes we have an idea on the maximal wave number 
{ r max ) which we want to resolve. For instance, if there is a source term in our equations with 
known band of frequencies. In this case we will solve (3.20) for a. Assuming that a is 
close to 1 we can simplify (3.20) and use instead 



. ( 'KT max 

\ 

N 

, ^moi ^ 2 ' 

(4.1) 


a = sm l N . 

) 

If 

_ N 
'T'ma.x 2 

j 

. N 
, j 

(4.2) 

then 

a = sin (| - 

f) 

- cos (n t ) ’ 

(4.3) 


j being the number of modes which we ‘give up’ resolving. Expanding in Taylor series 



which satisfies (2.19) and by using (4.4) and (2.20) we get 




2 1 
j + VP -FT N 


+ o (TV- 2 ) . 


(4.5) 


Remark 4-1 - Resolution analysis is closely related to maximal spacing analysis. By the 
sampling theorem we know that for any sampling interval A, there is maximal mode w c 
called the Nyquist critical frequency and is given by w c = y. For Fourier method we have 
Ax = ^ , hence the maximal mode which can be resolved is y. Similarly, in the Chebyshev 
case Ax max = ^ and the maximal mode is y which is equivalent to stating that 7r 
points per wavelength are needed for resolution. This result is given also in [5] based on 
expanding the trigonometric functions in Chebyshev polynomials. By Lemma 2.1 we get 
that for the transformed interpolating points, the maximal spacing is attained in the center 
of the interval. Therefore 


Ax 


max 


Substituting (4.3) in (4.6) we get 


S [ cos (l;^)]- J(0) 

arcsin[Qsin(j^)] 

arcsin(a) 


lim — — 

N— KX> /\x r , 


N 


-3 


as in (4.2). 


(4.6) 


(4.7) 


Accuracy considerations - For given e and TV, we can solve (3.36) for a and get an 
explicit expression 

«=7 TTT ’ (4.8) 


t + t~ l 

To examine the minimal spacing dictated by this choice, we expand a in Taylor series 

2 


Q = 1 - 1 ln! w (^) + 0 ( W_J ) • 


The expansion (4.9) is of the form (2.19), using it and (2.20) we get 

2tt 1 


Ax m i n = 


(|ln(e)| + \A- 2 +ln 2 (e)) 


N 


+ 0{N~ 2 ). 


(4.9) 


(4.10) 


The agreement between £ and the computed accuracy will depend on the constant 
c. Observe that there is a ‘give and take’ relation between resolution and accuracy. By 
decreasing r mox ,a is decreased (4.1) and therefore e is getting smaller (3.36). Hence, by 
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sacrificing the resolution of the high modes we get in return higher accuracy on the rest of 
the modes. See numerical results in Section 6. 


Adaptive approach - We have described above two formulas which give explicit expres- 
sions for a. One can also consider a third approach, an adaptive algorithm for computing 
a. Observing (2.6) we can regard the method described in this paper as a preconditioning 
one. For a given function with large gradients, we are looking for a parameter a such that 
after the transformation f(y) will be a smooth function. One can consider the tail of the 


series of pseudospectral Chebyshev coefficients to measure the smoothness of f(y)- Since, by 
stability considerations we would like a to be as large as possible, the adaptive algorithm 
should find N 

f | uU-k k( Q ) l 

otmax — max \ CL j / \i 

l £,=o l a «( a )l 

where e 0 is a given tolerance and a, are the computed Chebyshev coefficients. When the 
adaptive approach is implemented in time dependent problems, the search for an optimal 
a should restart whenever the solution behavior has been changed significantly. 


< e 0 , k < N 


} 


(4.11) 


x=g(y, <*,/))=- (g(y, a, p)-b), 


a 


5. Non-Symmetric Transformation 

The transformation function (2.3) is symmetric. The interpolating points (2.4) are dis- 
tributed symmetrically around the origin. When there is a boundary layer on one side of 
the domain, we would like to have the flexibility of putting more points on this side. To this 
end we modify the transformation (2.3) and take 

(5.1) 

(5.2) 

(5.3) 

(5.4) 

(5.5) 


where 


and 


g(y, a, /?) = arcsin 


1 


( 2at(fy + a — P 
V « + /? 


: ) 


b = - {<7(1; a,/?) + g {~ !;«»/?)} • 


For the derivative we have 


?(y; <*>/?) v«p 


and the parameters a and /? control the distribution of interpolation points near x - 1 
and x = — 1 respectively. An elaborate discussion of this transformation will be given in a 

future paper. 
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6. Numerical Results 


The following notations are used in this section: 

N = number of interpolating points, 

E = relative error of the derivative in the maximum norm, 
j = number of unresolved modes (4.2), 
yi = Chebyshev interpolating points (2.5), 

Z{ = check points, 

*i = 9 (cos(y, + i;a)) , 0 < i < N - 1. (6.1) 

In the first table we present the spectral radius p of AD (2*7) where cx. is given by (4,3). 
The spectral radius of the Chebyshev pseudospectral differentiation operator D is given in 
the last column. Using the new method for time dependent problems we have observed that 

A t(new method ) p{D) 

A t{Ch^Ih^j * p{AD)' 

Thus, from Table I we see that for N = 128 for example, the time step restriction of 
the new algorithm is almost 8 times larger then the one used in a standard Chebyshev 
discretization. 


Table I Spectral Radius 


j 

N 

a 

p{AD) 

P(D) 

1 

i'6 

0.9808 

18.927 

23.560 

1 

32 

0.9952 

42.897 

91.559 | 

1 

64 

0.9987 

92.286 

363.779 

1 

128 

0.9997 

192.165 

1452.706 

2 

16 

0.9239 

17.934 

23.560 

2 

32 

0.9808 

41.061 

91.559 

2 

64 

0.9952 

89.957 

363.779 

2 

128 

0.9988 

189.454 

1452.706 

3 

16 

0.8312 

18.624 

23.560 

3 

32 

0.9569 

44.600 

91.559 

3 

64 

0.9892 

97.846 

363.779 

3 

128 

0.9973 

204.997 

1452.706 


The results given in Tables II and III demonstrate the resolution and accuracy properties 
of the new algorithm and clarify the give and take 5 relation between the two as mentioned 
in Section 4. We applied the new differentiation algorithm to the trigonometric functions 

u k {x) = cos(jhnr) 1 < k < 16, (6.3) 
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and obtained approximations v' k , k = 1, • • • , 16. a is given by (4.3). Corresponding results 
for Chebyshev method are shown in the last column. In the last row we have printed £ as 
given by (3.36). The resolution property of f -j modes (4.2), (4.7) is clearly demonstrated 
in Tables II, III. As j increases, so does the accuracy of the modes resolved. As shown in 
Section 3 the new method should be exact if 


m = 


2kN 
N - 2j ' 


(6.4) 


is an even integer. This explains the high accuracy exhibited in relevant entries of Tables II 
and III. Comparing Table III to Table II we see that the accuracy by which the low modes 
are resolved is almost the same. The effect of increasing N is in the number of modes 
resolved with accuracy imposed by the choice of j. 


Table U N = 32 


k 

E(j = 1) 

E(i = 2) 

II 

■<-> 

EU = 4) 

E c h b 

1 

9.35E-03 

7.14E-04 

3.78E-05 

1.57E-06 

1.57E-12 

2 

1.86E-02 

1.35E-03 

6.23E-05 

1.77E-06 

5.83E-13 

3 

2.84E-02 

1.82E-03 

5.84E-05 

1.66E-06 

3.29E-13 

4 

3.70E-02 

2.02E-03 

2.06E-05 

2.97E-06 

5.65E-11 

5 

4.76E-02 

1.88E-03 

5.14E-05 

4.42E-06 

4.41E-08 

6 

5.44E-02 

1.25E-03 

1.42E-04 

3.90E-14 

7.56E-06 

7 

6.29E-02 

3.19E-14 

2.10E-04 

1.36E-05 

4.14E-04 

8 

7.04E-02 

2.04E-03 

1.68E-04 

2.93E-05 

9.33E-03 

9 

7.76E-02 

5.07E-03 

1.54E-04 

4.10E-14 

9.34E-02 

10 

8.36E-02 

9.26E-03 

1.07E-03 

2.46E-04 

4.52E-01 

11 

8.58E-02 

1.48E-03 

3.12E-03 

1.16E-03 

8.94E-01 

12 

8.46E-02 

2.10E-02 

6.60E-03 

3.01E-14 

1.56E+00 

13 

7.82E-02 

2.48E-02 

3.52E-14 

3.75E-01 

1.72E+00 

14 

5.70E-02 

4.07E-14 

7.39E-01 

1.32E+00 

1.70E+00 

15 

3.87E-14 

1.25E+00 

1.79E+00 

1.86E+00 

1.34E+00 

16 

1.69E+00 

1.92E+00 

1.76E+00 

1.79E+00 

1.63E+00 


4.32E-02 

1.79E-03 

6.99E-05 

2.49E-06 

— 
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Table in N = 64 


k 

E{j = 1) 

II 

E(] = 3) 

E(j = 4) 

E c i ib 

1 

4.55E-03 

3.53E-04 

2.04E-05 

1.04E-06 

3.49E-11 

3 

1.36E-02 

1.03E-04 

5.53E-05 

2.49E-06 

7.16E-12 

5 

2.26E-02 

1.60E-03 

7.29E-05 

2.26E-06 

3.50E-12 

| 7 

3.16E-02 

1.99E-03 

6.35E-05 

1.99E-13 

1.85E-12 

9 

4.04E-02 

2.13E-03 

2.14E-05 

3.53E-06 

1.02E-12 

11 

4.90E-02 

1.93E-03 

5.26E-05 

6.13E-06 

2.15E-12 

13 

5.74E-02 

1.26E-03 

1.47E-04 

4.12E-06 

1.55E-08 

15 

6.52E-02 

6.58E-03 

2.32E-04 

6.77E-06 

1.48E-05 

17 

7.25E-02 

1.99E-03 

2.44E-04 

2.76E-05 

2.62E-03 

19 

7.89E-02 

4.89E-03 

6.26E-05 

4.56E-05 

9.76E-02 

21 

8.38E-02 

8.91E-03 

5.40E-04 

5.81E-14 

6.61E-01 

23 

8.66E-02 

1.42E-03 

1.96E-03 

3.12E-04 

1.78E+00 

25 

8.57E-02 

2.07E-03 

4.94E-03 

1.53E-03 

1.76E+00 

27 

7.85E-02 

2.72E-03 

1.01E-02 

4.48E-03 

1.86E+00 

29 

5.82E-02 

2.46E-02 

9.01E-14 

7.12E-01 

1.62E+00 

31 

1.48E-13 

1.51E+00 

2.03E+00 

1.91E+00 

1.77E+00 

ILx_ 

4.32E-02 

1.84E-03 

7.80E-05 

3.21E-06 

— 


In Tables IV and V we present mesh refinement charts for the functions 

pf \ 0-05 

^ ~ X 2 + 0.05 ’ 


(6.5) 


/*(*) = 


exp(2x) 


( 6 . 6 ) 


2 + cos(15x) 

respectively. In Table IV, a is given by (4.3) with j = 3. Observe the fast convergence up 
N 64 . The error is not decreasing beyond this points, since all the modes have been 


resolved to the accuracy enforced by the choice of j. In fable V, a is computed by (4.8) 
with £ l.E 05 . The results for Chebyshev method are given in the last column. 


Table IV Mesh refinement chart, fi(x) = x °±q 0S 


N 

a 

E(j = 3) 

E c h b 

16 

0.83147 

9.394E— 02 

1.777E— 01 

32 

0.95694 

1.019E-03 

9.281E— 03 

64 

0.98918 

1.507E— 06 

1.486E— 05 

128 

0.99729 

1.794E— 06 

8.845E— 11 

256 

0.99932 

1.939E— 06 

9.923E— 11 
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Table V Mesh refinement chart, .^(x) = 2 +cL\ 2 i 5 j;) 


N 

a 

E(e = l.£-05) 

&chb 

16 

0.80761 

1.853E— 01 

1.717E— 01 

32 

0.94208 

9.007E— 02 

6.461E— 02 

64 

0.98452 

2.298E— 03 

9.032E— 03 

128 

0.99603 

1.968E— 06 

5.774E— 05 

256 

0.99899 

6.344E— 06 

1.669E— 09 


We have solved the model problem described in the introduction (1.1)- (1.3) and the 
results are reported in Table VI. The solution is computed at t — 1 using fourth order 
Runge-Kutta as time marching algorithm. The initial and boundary conditions are 


u°(x) = [xexp cos (m7rx) - (— l) m ] , 


s(t) = 0 


(6.7) 

(6.8) 


respectively. The numerical solution is compared to the exact solution 


u(x,t) = 

The results presented in the table provide a comparison between the new method, where a 
is computed by (4.3) with j = 1, and standard Chebyshev method, nsteps is the number 
of time steps. For m = 6 we needed 91 points, in the Chebyshev case, in order to 
achieve the accuracy given in the last column. For stability we had to use 1300 time steps . 
Taking smaller At would not decrease the error as shown in the second row. Using the new 
algorithm, we took only 65 points. 80 time steps were sufficient for stability. In order to 
get accuracy close to the one we had in the Chebyshev case, 200 time steps were needed. 
Reducing At beyond this point would not reduce the error significantly as shown by the 
fifth row. In the next set of experiments we took m = 12. In both cases, we had to double 
the number of points in order to resolve the solution. The results are provided in the rest of 
the table. 


x + t > 1 


[{x + t ) exp-^ + ‘- 1 > 2 cos[m7r (x + t )] - (-1)” 1 ] x + t < 1 . 


(6.9) 
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Table VI Time-dependent problem, (l.l)-(l.S) 


m 

N 

a 

nsteps 

max(|i£|)/ max(|ii|) 

6 

91 

0 

1300 

8.59E-03 

6 

91 

0 

2600 

8.53E-03 

6 

65 

0.9988 

80 

1.79E-01 

6 

65 

0.9988 

200 

7.08E-03 

6 

65 

0.9988 

300 

5.64E-03 

12 

181 

0 

5200 

5.32E-03 

12 

129 

0.9997 

200 

2.22E-01 

12 

129 

0.9997 

600 

5.71E-03 


Lamb Problem 

The problem is of wave propagation in a uniform and isotropic elastic two dimensional 
halfspace subjected to a point source applied in the vicinity of the free surface. This problem 
is numerically challenging because of. the presence of Rayleigh surface waves around the 
free surface the calculation of which requires an accurate representation of the boundary 
conditions. 

Let x and y denote horizontal and vertical cartesian coordinates respectively and t 
the time variable. The system of equations to be solved is 


dV x 

dt 

1 / d<7 X X d(T X y \ 

P\dx dy ) +Jx 

(6.9) 

dV„ 

dt 


(6.10) 

dcr xx 

dt 

- 

(6.11) 

dOyy 

dt 

,dV x n , dV y 

= X 3 X ^ X + 2 ^ay 

(6.12) 

der xy 

dt 

dV, , t>V„ 
= “dy +>, dz- 

(6.13) 


V x and V y respectively denote the horizontal and vertical velocities, a xx , a yy and a xy 
are the stress components, f x and f y are the body forces, p is the density and A and 
p are Lamb’s constants. The system is the same as the one used by Bayliss et. al. [1] for a 
fourth order finite difference scheme. 

During the calculations, the variables 14, V y , cr xx , a yy and a xy are advanced in time 
after specification of the body forces f x and f y . In this work we choose to approximate the 
horizontal derivative by the Fourier method whereas for the vertical coordinate y we choose 
the modified Chebyshev discretization as described in this paper, using the transformation 
function (5.1). The boundary conditions at y = 0 is cryy = <r xy = 0 whereas for the 
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bottom boundary y = L we choose the condition that the incoming characteristics are 
zero [1]. In addition an absorbing strip was applied along the lower boundary and the sides 
of the grid to prevent reflections or wraparound from the boundaries [7], For the present 
problem the material parameters had uniform values of p = 1 •2 jr / cm 2, and P and S 
velocities of V p = = 2000 m/sec and V 3 = ^ = 1155 m/aec respectively. For the 

body forces, f x = 0 and f y = £(x — xo)$(y — yo)h(t ), where xo = 250 m , yo — l-8 m 
and h(t) was a band limited Ricker wavelet with highest frequency of 40Hz [11]. For 
the spatial discretization Ax = A y = 10 m and the grid modification parameters (5.1) 
are a = 0.729 and /? = 0.620. The solution was advanced in time to l aec hy the 
fourth-order Runge-Kutta method using a time step of 0.002 seconds. This time step is 
approximately seven times larger then the maximum allowable time step for an ordinary 
Chebyshev discretization, and is approximately equal to the time step which would be used 
with a uniform Fourier grid from accuracy considerations (e.g. ^ w 0.4 ). 

Figures 1—3 present a comparison between the numerical and analytical horizontal dis- 
placement time histories at points located at respective distances of 200 m ,500 m and 800 m 
from the source. A corresponding comparison of vertical displacements is presented in Fig- 
ures 4—6. As the figures show the match between numerical and analytical solutions is 
virtually perfect. 
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APPENDIX 


Lemma A.l: For m even 


cos (mVO = (-l)“r m (sin(V>)). 

For m odd 

sin(mV’) = (— l) c ^ J 'T m (sin(V’)). 


(A.l) 
(A- 2) 


Proof: Chebyshev polynomials satisfy the recurrence relation 

r n+ i(s) = 2xT n (x) - T n -i(x). (A. 3) 

Using basic trigonometric identities, the reccurence relation and induction we get 

cos((m + 2)V>) = 2[2 cos 2 (^>) — 1] cos (mtp) — cos[(m — 2)ip] 

= (-1) “ +1 {2[2sin 2 (^) - l]T m (sin(^)) - T m _ 2 ( sin(V>))} 

= (-l)“ +1 T’ m+2 (sin(^)). (A. 4) 

We will use now (A. 4) to show (A. 2). 


sin((m + 2)xp) 


cos(2i/>) sin(m^) + cos(m^) sin(2^) 

(l — 2 sin 2 sin(^)) + sin(^>) {cos ((m + 1)^>) + cos ((m — 1)^)} 

(-l) 21 ^ j(2 sin 2 xf) — l)T ro ( sin(V’)) + sin(^) [T m+1 (sin(^)) - T m _ 1 (sin(V’))]} 
(-l) 22 ^ {2 (2 sin 2 (i/>) - l) r m (sin(V>)) - T m _ 2 (sin(i/>))} 

(-l)^T m+2 . (A. 5) 


Lemma A. 2: For m odd 


cos (imp) — cos(^)P TO _i(sin(V»)) 


(A.6) 


and for m even 


sin(mi/>) = cos(V>)Qm-i(sin(VO) (A.7) 

where P m -\, Q m - 1 are polynomials of degree m — 1. 

Lemma A. 2 is easily verified by using trigonometric identities, induction and the results of 
Lemma A.l. 
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Figure 1. Comparison between the numerical and analytical horizontal displacement 

time history at distance of 200 m from the source. 


Leaend 



Comparison between the numerical and analytical horizontaJ displacement 
time history at distance of 500 m from the source. 




23 


Figure 3. Comparison between the numerical and analytical horizontal displacement 

time history at distance of 800 m from the source. 
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Figure 4. Comparison between the numerical and analytical vertical displacement 

time history at distance of 200 m from the source. 
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Legend 

Analytical solution 
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Comparison between the numerical and analytical vertical displacement 
time history at distance of 800 m from the source. 
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